Turing patterns in network-organized activator-inhibitor systems 

Hiroya Nakao* 

Department of Physics, Kyoto University, Kyoto 606-8502, Japan 

Phone: +81-75-753-3742 FAX: +81-75-753-3819 

E-mail: nakao@ton.scphys.kyoto-u.ac.jp 

O"' 

r— I , 

^^ , Alexander S. Mikhailov 

Department of Physical Chemistry, 

a ■ 

^H ■ Fritz Haber Institute of the Max Planck Society, 

pg . Faradayweg 4-6, 14195 Berlin, Germany 

Phone: +49-30-8413-5122 FAX: +49-30-8413-5106 

Q ■ E-mail: mikhailov@fhi-berlin.mpg.de 

(Dated: May 13, 2010) 

f^h \ Abstract 

Turing instability in activator-inhibitor systems provides a paradigm of nonequilibrium pattern formation; 

it has been extensively investigated for biological and chemical processes. Turing pattern formation should 

>' 
\Q ' furthermore be possible in network-organized systems, such as cellular networks in morphogenesis and 

00 '. 

ecological metapopulations with dispersal connections between habitats, but investigations have so far been 

restricted to regular lattices and small networks. Here we report the first systematic investigation of Turing 

patterns in large random networks, which reveals their striking difference from the known classical behavior. 

In such networks, Turing instability leads to spontaneous differentiation of the network nodes into activator- 

^. , rich and activator-low groups, but ordered periodic structures never develop. Only a subset of nodes having 

close degrees (numbers of links) undergoes differentiation, with its characteristic degree obeying a simple 

general law. Strong nonlinear restructuring process leads to multiple coexisting states and hysteresis effects. 

The final stationary patterns can be well understood in the framework of the mean-field approximation for 

network dynamics. 



I. INTRODUCTION 

In 1952, A. Turing published a seminal paper [l|, demonstrating that differences in diffusion 
rates of reacting species can alone destabilize the uniform state of the system and lead to the 
formation of spatial patterns and suggesting this as a possible mechanism of biological morphogen- 
esis. The Turing instability and resulting patterns have subsequently been theoreticallyanalyzed 
and experimentally confirmed for various chemical, biological, and ecological systems [2h6]; they 
are generally viewed as a paradigm of nonequilibrium pattern formation. In 1971, Othmer and 
Scriven [7j] pointed out that the Turing instability is also possible in network-organized systems and 
this should be important for the understanding of multi-cellular morphogenesis. At an early stage 
of the organism development, a network of inter-cellular connections is formed and morphogens 
are diffusively transported over such a network; differentiation of cells may thus be induced by 
the network Turing instability. On the other hand, many ecosystems represent metapopulations 
distributed over discrete habitats forming networks with dispersal connections [a, |9j. Prey and 
predator or host and parasite species may migrate over such networks with different diffusional 
mobilities. Recently, diffusional spreading of infectious diseases over airline transportation net- 
works has attracted much attention [10]. Epidemic dynamics on the networks can be described 
by reaction-diffusion models where infected and susceptible individuals play the role of interacting 
species [111 . [l2j. In all such network-organized systems with diffusional transport of interacting 
species, the Turing instabilities and resulting nonlinear patterns can generally be expected. 

Although nonlinear dynamics on complex networks is attracting significant attention, most 
of the investigations have focused on synchronization phenomena of oscillator networks (see, e.g., 
recent reviews [13|, [lj] ) . Despite the potential importance of network Turing patterns and the large 
amount of research on classical Turing patterns in spatially extended systems, very little research 
on this problem has been performed so far. Early studies by Othmer and Scriven [7|, [15] have 
provided abstract mathematical framework for the analysis of network Turing instability, but they 
explicitly considered only simple examples of regular lattices close in their properties to continuous 
media. Recently, Horsthemke et. al. [la . Il7| ] have discussed the possibility of Turing instability in 
coupled chemical reactors, but only for small networks. 

In this article, we present the systematic analytical and numerical study of the Turing instability 
and of the developing nonlinear patterns in large random networks. We find that the Turing 
instability generally occurs in network-organized activator-inhibitor systems, but its properties are 
very different from those characteristic for the classical continuous media. 



II. NETWORK-ORGANIZED ACTIVATOR-INHIBITOR SYSTEMS 



Classical activator- inhibitor systems in continuous media are described by equations of the form 

-u(x,t) = f(u,v) + D act V u(x,t), 



dt 

d 

—v(x,t) = g(u,v) + D inh V 2 v(x,t), (1) 

where u(x, t) and v (x, t) are local densities of the activator and inhibitor species. Here, the 
functions f(u,v) and g(u,v) specify dynamics of the activator that autocatalytically enhances 
its own production and of the inhibitor that suppresses the activator growth. D ac t and D^^ are 
the diffusion constants of the activator and inhibitor species. The classical Turing instability [l( 
sets in as the ratio Di n h/D act of the two diffusion constants is increased and exceeds a threshold. 
It leads to spontaneous development of alternating activator-rich and activator-low domains from 
the uniform background. 

In our study, we consider the network analog of the model ([I]), where activator and inhibitor 
species occupy discrete nodes of a network and are diffusively transported over the links connecting 
them. We consider a connected network consisting of N nodes, i = 1, ■ • • , N. The network topology 
is defined by a symmetric adjacency matrix whose elements A^ take values Aij = 1 if the nodes i 
and j are connected (i ^ j) and Aij = otherwise. Diffusive transport of the species into a certain 
node i is simply given by the sum of incoming fluxes to the node i from other connected nodes {j}, 
where the fluxes are proportional to the concentration difference between the nodes (Fick's law). 
By introducing the network Laplacian matrix whose elements are given by Lij = A^ — kiSij, where 
k{ = J2j=i Aij is the degree of the node i, diffusive flux of the species u to node i is expressed as 
J2j=i LijUj = J2j=i A{j(uj — U{), and similarly for v. Generally, diffusional mobilities of species u 
and v on the network are different. The equations describing network-organized activator- inhibitor 
systems are thus given by 

d N 

-jliH(t) = f(ui,Vi) + e'J2LijUj, 

d N 

-^7^(t) = g(ui,Vi) + ae^2LijVj, (2) 

i=i 
for % = 1, • • • , N. Here, f(u,v) and g(u,v) specify the local activator-inhibitor dynamics on indi- 
vidual nodes. We denote the diffusional mobility of the activator species as e(= D ac t) and of the 
inhibitor species as ae(= Di n h), so that a = Di n h/D ac t is the ratio between them. The considered 



systems have a uniform stationary state (u,v), where f(u,v) = and g(u,v) = 0. This uniform 
state can become unstable as a result of the Turing instability. If u and v correspond to the acti- 
vator and inhibitor species, functions / and g should satisfy several conditions which are given in 
the Methods section. 

As the examples of activator-inhibitor systems, we use the Mimura-Murray model of prey- 
predator populations [3] and the classical Brusselator model [2] which are described in the Methods 
section. This study is focused on the Turing instability and pattern formation in large random 

iaU,U,U 



networks. We use scale- free networ 
Erdos-Renyi random networks [18 



cs which are ubiquitous in Nature 

n 



and the classical 



i, llSI], both described in the Methods section. For convenience, 
network nodes {i} are always sorted below in the decreasing order of their degrees {ki}, so that 
the condition k\ > hi > • • • Ajjv holds. 

III. THE TURING INSTABILITY 

The Turing instability is revealed through the linear stability analysis of the uniform stationary 
state with respect to nonuniform perturbations (see Methods for the details). In the classical case of 
continuous media [l|, nonuniform perturbations are decomposed over a set of spatial Fourier modes, 
representing plane waves with different wavenumbers. As has been originally noticed by Othmer 
and Scriven 7|], in the networks, the role of plane waves should be played by eigenvectors of their 
Laplacian matrices. The ei gen values A a and eigenvectors <p Q ) = (<f)\ , • • • , <Pn ) °^ ^ ne Laplacian 
matrix Ly are determined [l3|, l20(] by Y^j=i ^ij4>j = h~ a (j)\ , with a = 1, • • • ,N. All eigenvalues 
of Lij are non-positive. We sort the indices {a} in the decreasing order of the eigenvalues, so that 
the condition = Ai > A2 > • • • > A^v always holds. 

Introducing small perturbations (Sui, Svi) to the uniform state as (ui,Vi) = (u,v) + (5v,i,5vi) 
and substituting this into equations ([2]), a set of coupled linearized differential equations is ob- 
tained. By expanding the perturbations over the set of Laplacian eigenvectors as 5ui(t) = 
J2a=i c a ex P [A Q £] (t>i an d Svi(t) = J2a=i c aB a exp [X a t] <j>f , these equations are transformed into 
N independent linear equations for different normal modes. The linear growth rate X a of the a-th 
mode is determined from the characteristic equation. When Re X a is positive, the a-th mode is 
unstable. The Turing instability takes place when one of the modes (i.e., the critical mode) begins 
to grow. At the instability threshold, Re X a = for some a = a c and Re X a < for all other 
modes. In the Turing instability, the critical mode is not oscillatory, Im A Qc = 0. 

As an example, Fig. [TJshows the growth rate A as a function of A for the Mimura-Murray model. 



At e = 0.060, three curves corresponding to different ratios a of diffusion constants (below, at and 
above the instability threshold) are displayed. In this figure, critical curves for two other values of 
the parameter e are also shown. 

Generally, the Turing instability becomes possible for a > a c . The dispersion curve A = 
-F(eA) first touches the horizontal axis at A = A c and the Laplacian mode (j)\ a c , possessing the 
Laplacian eigenvalue A Qc that is closest to A c , becomes critical. For the critical mode, the coefficient 
B a is positive, so that when the activator concentration increases, the inhibitor concentration 
also increases accordingly. Explicit expressions are given in the Methods section. Note that the 
Laplacian spectrum of a network is discrete and, therefore, the instability actually occurs only 
when one of the respective points on the dispersion curve crosses the horizontal axis. 

The above results are analogous to those holding for continuous media (cf. 2l|). The critical 



ratio <7 C in the networks is the same as in the classical case. The Laplacian eigenvalue A c of the 
critical network mode corresponds to — q%, where q c is the wavenumber of the critical mode in the 
continuous media. Despite such formal analogies, properties of the network Turing patterns are 
very much different from their classical counterparts, as demonstrated in the following sections. 

IV. LOCALIZATION OF LAPLACIAN EIGENVECTORS AND CHARACTERISTIC 
PROPERTIES OF CRITICAL TURING MODES 

When a Turing pattern starts to grow after slightly exceeding the instability threshold, the 
activator distribution in this pattern is determined by the critical Laplacian eigenvector, i.e. we 
have 5ui oc % • Therefore, to understand organization of the growing Turing patterns, properties 
of Laplacian eigenvectors should be considered. 

As an example, Figs. [2ja,b) display critical eigenvectors of a network for two different values of 
the diffusion constant e. The same eigenvectors are shown graphically in Figs.[2jc,d). In the chosen 
representation, network nodes with larger degrees (hubs) are located in the center and the nodes 
with lower degrees in the periphery of the graph. The nodes are colored red when (pf > 0.1 (e.g. 
activator concentration is significantly increased), blue when (pf* < —0.1 (significantly decreased), 
and yellow for —0.1 < 4>i < 0.1 (no significant change). 

It is clearly seen in Fig. [2] that spontaneous differentiation of the nodes takes place - the dis- 
tinguishing feature of the Turing instability. However, it affects only a fraction of all nodes. The 
differentiated nodes, with significant deviations of the activation level, tend to have close degrees. 
When the diffusional mobility e is small, only a subset of hub nodes undergoes differentiation 



[Figs.[2ja),(c)]. If e is large, differentiated nodes have just a few links [Figs. 12(b), (d)]. Thus, a cor- 
relation between the characteristic degrees of the differentiated nodes and the diffusional mobility 
is present. The behavior observed in Fig. [2jis general. As we show below, it is related to the effect 
of localization of Laplacian eigenvectors in large random networks. 
As has recently been discussed 



22 ] . Laplacian eigenvectors in large random networks with 
relatively broad degree distributions tend to localize on subsets of nodes with close degrees. The 
localization effect for a scale-free network is illustrated in Fig. [3j Here, all nodes are divided into 
groups with equal degrees k. For each group k and a given Laplacian eigenvalue A, the number 
of "differentiated" nodes with cpf 1 > 0.1 or cpf 1 < —0.1 in the respective eigenvector (pf is 
counted. The density diagrams in Fig. 3 display in the color code the relative numbers of such 
nodes as functions of the Laplacian eigenvalue A and the degree k. Examining Fig. [3j one can see 
that differentiated nodes are approximately located along the diagonal of the density map. The 
localization effect is more pronounced for the larger network of size N = 1000. Similar localization 
behavior is observed for the Erdos-Renyi networks (see Supplementary information). Thus, we 

(a) 

see that each Laplacian eigenvector (j)\ is characterized by its characteristic localization degree 
k a . Moreover, this characteristic degree is approximately equal to the negative of the respective 
eigenvalue, so that a simple relationship k a ~ — A Q holds. 

On the other hand, as implied by the linear stability analysis (see Methods), the growth rate X a 
of each mode depends only on the combination eA a of the diffusional mobility e and the eigenvalue 
A a of that mode, i.e. we have X a = F(eA a ). Therefore, the Laplacian eigenvalue A Qc of the critical 
mode a c with A Qc = should be inversely proportional to the diffusive mobility e, i.e., A ac ex 1/e. 
Hence, modes with large negative eigenvalues A a tend to become unstable for the small mobilities 
e (note that A a < in our definition). 

Combining the two relationships, k a ~ — A Q and A ac oc 1/e, a simple scaling law k ar oc 1/e 
is obtained. It implies that the characteristic degree k ac of the differentiating subset is inversely 
proportional to the diffusional mobility e. 

The dependence A Qc oc 1/e holds for any activator-inhibitor model exhibiting the Turing in- 
stability. The localization of Laplacian eigenvectors has been observed by us (to be separately 
published) also for other random networks with relatively broad degree distributions. The charac- 
teristic localization degree k ac of the critical Turing mode is generally a monotonously increasing 
function of the negative of the critical Laplacian eigenvalue, — A Qc , and thus a decreasing function 
of the diffusional mobility e. 



V. STATIONARY TURING PATTERNS 

The initial exponential growth is followed by a nonlinear process, leading to the formation of 
stationary Turing patterns. The nonlinear evolution of the system and the properties of emerging 
stationary patterns have been studied by us in numerical simulations. Figure H] presents typical 
results, obtained for intermediate diffusional mobility (e = 0.12) and slightly above the instability 
threshold (a = 15.6) for the Mimura-Murray model on the random scale-free network of size 
N = 1000 and mean degree (k) = 20. The nodes are sorted in the order of their degrees, as shown 
in Fig. [^d). 

Starting from almost uniform initial conditions with small perturbations, exponential growth is 
observed at the early stage. The activator pattern at this stage, Fig. Sl^b), is similar to the critical 
mode, Fig. [U^a), with the deviations due to the contributions from neighboring modes that are 
already excited to some extent. Later on, however, strong nonlinear effects develop, and the final 
stationary pattern, Fig. HJc), becomes very different from the one determined by the critical mode. 

Observing the nonlinear development, we notice that some nodes get progressively kicked off the 
main group near the destabilized uniform solution in this process (see Video in the Supplementary 
information). Eventually, in the asymptotic stationary state, the nodes become separated into two 
groups. The separation into two groups occurs only for the nodes with relatively small degrees 
(roughly i > 200, k{ < 24), while the nodes with high degrees {i < 200, k{ > 24) do not undergo 
the differentiation. 

Our numerical investigations furthermore reveal that the outcome of nonlinear evolution de- 
pends sensitively on initial conditions. Different Turing patterns are possible at the same parameter 

values and strong hysteresis effects are observed. As an example, Fig. [5ja) shows how the amplitude 

r 1 1/ 2 

of the stationary Turing pattern, defined as A = X)i=i {( u i ~ ^) 2 + ( v i ~ ^) 2 } : varies under 

gradual variation of the parameter a in the upward or downward directions. Stationary patterns 

observed at points P, Q, and R in Fig. Ufa) are shown in Fig. \Mjo)- 

As a was increased starting from the uniform initial condition, the Turing instability took place 

at a = a c , with the amplitude A suddenly jumping up to a high value that corresponds to the 

appearance of a kicked-off group. If a was further increased, the amplitude A grew. Starting to 

decrease a, we did not however observe a drop down at a = a c . Instead, a punctuated decrease in 

the amplitude A, which is characterized by many relatively small steps, was found. Reversing the 

direction of change of the parameter a at different points, many coexisting solution branches could 

be identified. The characteristics of Turing patterns vary with their amplitudes. When A is close 



to zero [ point R in Fig. [5ja) ], only a few kicked-off nodes remain in the system. Such localized 
Turing patterns with only a small number of destabilized nodes can coexist with the linearly stable 
uniform state and are found below the Turing instability threshold, for a < o~ c . 

VI. THE MEAN-FIELD THEORY 

To understand properties of the developed Turing patterns above the instability boundary 
(o~ > <7 C ), one can use the mean-field approximation, similar to that previously employed in the 

n n 

investigations of epidemics spreading on networks [23] and for the networks of phase oscillators [24] . 
We start by writing equations (|2|) in the form 

— Ui{t) = f(ui,Vi) +e(/i- -kiUi), 

-rVi(t) = g(ui,Vi) + ae(hl v ' -hvi), (3) 

where local fields felt by each node, h\ u = J2j=i^ij u j an d nf = J2j=iAijVj, are introduced. 
These local fields are further approximated as h\ — kiH^ u > and h\ — kiH^ v \ where global mean 
fields are defined by H^ u > = Ylj=i w j u j an d H^ v > = J2j=i w j v j- The weights Wj = kj/ [J2e=i kt 
ki/ktotai take into account the difference in contributions of different nodes to the global mean 
field, depending on their degrees (cf. [23j, |24|). Thus, the local fields are taken to be proportional 
to the degree of a node, ignoring the details of its actual connections. 

With this approximation, the individual activator-inhibitor system on each node interacts only 
with the global mean fields, and its dynamics is described by 

j t u{t) = f(u,v) + P(HM-u), 

j t v{t) = g(u,v)+o-(3(H^-v). (4) 

We have dropped here the index i, since all nodes obey the same equations, and introduced the 
parameter j3{i) = eki. If diffusion ratio a is fixed and the global mean fields H^ u ' and H^"' are 
given, the parameter /3 plays the role of a bifurcation parameter that controls the dynamics of 
each node. Equations @ have a single stable fixed point when j3 = (i.e. e = 0), and, as f3 is 
increased, this system typically undergoes a saddle-node bifurcation that gives rise to a new stable 
fixed point. 

As an example, we have computed stationary Turing patterns for the Mimura-Murray model by 
numerical integration of equations ([2]) and determined the respective global mean fields H^ u > and 



H( v > at a = 15.6 and a = 30. Substituting these computed global mean fields into equations (|3J), 
bifurcation diagrams for a single node have been obtained (solid curves in Fig. [B] (a,c)). In this 
example, one of the two stable branches vanishes by another saddle-node bifurcation when (3 is 
increased further. These diagrams can be compared with the actual stationary Turing patterns. 
Each node i in the network is characterized by its degree k%, so that it possesses a certain value 
of the bifurcation parameter, (3 = ski. Therefore, the Turing pattern can be projected onto these 
bifurcation diagrams, as shown by crosses in Fig. [6|a,c). We see a relatively good agreement 
between the stable branches and the data from the actual Turing patterns. Furthermore, we 
directly compare in Fig. [6^b,d) the computed Turing patterns with the mean-field predictions, 
based on equations §%§. The Turing patterns are well fitted by the stable branches, though the 
scattering of numerical data gets enhanced near the branching points. 

In the Supplementary information, a similar analysis is performed for the Brusselator model. 
This model has a different bifurcation diagram in the presence of external fields. Nonetheless, 
a good agreement with the predictions of the mean-field theory is again observed. Thus, fully 
developed Turing patterns in the networks are essentially explained by the bifurcation diagrams of 
a single node coupled to global mean fields, with the coupling strength determined by the degree 
of the respective network node. The mean-field theory is generally not applicable for the localized 
Turing patterns found below the Turing instability threshold. 

VII. DISCUSSION AND CONCLUSIONS 

The fingerprint property of the classical Turing instability in continuous media is the sponta- 
neous formation of periodic patterns above the critical point. Our investigations of the Turing 
problem for large random networks have revealed that, while the bifurcation remains essentially 
the same, properties of the emerging patterns are strongly different. In the networks, the critical 
Turing mode is localized on a subset of network nodes with the degrees close to a characteristic 
degree controlled by the mobility of species. The final stationary patterns are much different from 
the critical mode. Multistability, i.e. coexistence of a number of various stationary patterns for 
the same parameter values, is typically found and the hysteresis phenomena are observed. Above 
the instability threshold, network Turing patterns can be well understood in the framework of 
the mean-field approximation. In this approximation, each network element is coupled to certain 
global fields collectively determined by the entire system, and interactions of the element with its 
neighbor nodes are neglected. The strength of coupling to the global fields depends on the number 



of links connecting an element to the rest of the network. 

The Turing instability is also possible in regular lattices representing a special case of networks. 
An activator-inhibitor system on a lattice can be viewed as a finite-difference approximation for 
the respective reaction-diffusion problem in the space and Turing patterns on the lattices should 
therefore have almost the same properties as in the continuous media. The divergent behavior, 
characteristic for large random networks, must be related to a difference in the structural properties 
of such systems. Diameters of random networks are relatively small and nodes in such networks 
cannot be separated by large distances (roughly estimated as L = IniV for the Erdos-Renyi and 



scale- free networks 25f] ). For comparison, a cubic lattice with N nodes in the d-dimensional space 
has the diameter about L = N 1 ^. Thus, a lattice with the same size iVasa random network and 
a comparable diameter L must have a high dimension d S> 1. In lattices with high dimensionality 
and short lengths, Turing patterns with many alternating domains are however not possible and 
just a few domains (clusters) shall be found, resembling what is indeed observed by us in large 
random networks. 

Because of their small diameters, diffusional mixing in random networks is strong. Large random 
networks are, therefore, structurally much closer to the globally coupled systems than to the 
low-dimensional lattices. Globally coupled activator-inhibitor populations have previously been 
considered and spontaneous separation of the elements into two groups has also been found in such 
systems 



261 ] . There is, however, an important further aspect distinguishing random networks from 
the lattices or simple globally coupled systems. All nodes in a lattice (or in a globally coupled 
population) are equivalent and have the same degree (number of neighbors). In contrast to this, 
random networks effectively represent strongly heterogeneous systems. They are characterized 
by broad degree distributions (less broad but still relatively wide for the finite-size Erdos-Renyi 
networks). This heterogeneity becomes essential in the problems involving diffusion. Under the 
same concentration gradients across the links, a node with a higher degree receives a larger incoming 
flux from the neighboring nodes. 

The heterogeneity is responsible for the localization of Laplacian eigenvectors on the subsets 
of nodes with close degrees. Laplacian eigenvectors of networks are known to play an important 
role in the synchronization phe nomena. However, only the second and the last of such eigenvectors 



are significant there (see [13|])- In contrast, in network Turing problems, the entire Laplacian 
spectrum becomes significant. By varying the diffusional mobility of species, critical Turing modes 
corresponding to different Laplacian eigenvectors are realized. 

In the present study, a general framework for the analysis of network Turing patterns has been 

10 



proposed. Numerical investigations, confirming the theory, have been performed for the ecological 
predator-prey Mimura-Murray model [3J and for the classical chemical Brusselator model [2fl • Both 
models belong to the activator-inhibitor class. There is moreover a different broad class of models 
where the first species is characterized by autocatalytic growth (i.e., represents an activator) and 
it consumes for its growth the second species which effectively represents a renewable resource 
(see, e.g., [2l|]). In these systems, growth of the activator leads to the depletion of the renewable 
resource, which has an inhibitory effect on the activator. The considered Turing instabilities should 
also exist in such other network-organized systems. 

The results of our study may be important in a broad range of applications. Turing instabilities 
can generally be expected in various cellular, ecological or epidemic networks in Nature and their 
detection and observation represent a major scientific challenge. With the progress in engineering of 
synthetic ecosystems 



27[ , artificial ecological networks exhibiting Turing patterns can be designed 



in the future. 



VIII. METHODS 



Since u is the activator and v is the inhibitor in equations ([2]), partial derivatives of f(u, v) and 
g{u,v) at (u,v) should satisfy the following conditions: f u = df/du\, Ui} \ > 0, f v = df/du\, u ^ < 0, 
9u = 9g/du\t us \ > 0, and g v = dg/dv\i^^. < 0. The uniform stationary state of the system 
(ui,Vi) = (u,v) for all i = 1, • • • ,N is assumed to be linearly stable in the absence of diffusion, 
which requires /„ + g v < and f u g v - f v g u > 0. 

In the Mimura-Murray model [3j], u and v correspond to the prey and the predator densities. 
In this model, we have f(u,v) = {(a + bu — u 2 )/c — v}u and g(u,v) = {u — (1 + dv)}v, where 
the parameters have been chosen as a = 35, b = 16, c = 9, and d = 2/5 in the present study, 
yielding the fixed point (u,v) = (5, 10). In the Brusselator model [2], variables u and v correspond 
to densities of chemical activator and inhibitor species. Here we have f(u, v) = p — (r + l)u + u 2 v 
and g(u, v) = ru — u 2 v and the parameters have been chosen as p = 1 and r = 1.8. The fixed point 
is (u,v) = (1,1.8). 

Scale-free networks are generated by the preferential attachment algorithm of Barabasi and 
Albert [la. Il9j]. in which nodes with larger degrees tend to acquire more links. Starting from m 
fully connected initial nodes, we are adding m new connections at each iteration step, so that the 
mean degree is (k) — 2m. The simple Erdos-Renyi networks are generated by preparing TV" nodes 
and then randomly connecting two arbitrary nodes with probability q, yielding the mean degree of 



11 



(k) ~ Nq [19]. 

The Laplacian Ljj of any network is a real, symmetric, and negative semi-definite matrix [20] . All 
eigenvalues are real and non-positive. The eigenvectors are orthonormalized as J2i=i 4>i 4% = ^a,/3 
where a, /3 = 1, • • • , N. 

The linear stability analysis is performed in close analogy to the classical case of continuous 
media. Small perturbations Sut and Svi obey linearized differential equations 

N 

dt 



-5ui(t) = f u 5ui + f v Svi + e ^2 LijSiij, 



d N 

—Svi(t) = g u 5ui + g v 5vi + ae ^ Lijdvj. (5) 

Expanding Sui and 5vi over the Laplacian normal modes (pf as described in the main text, the 
following eigenvalue equation is obtained: 

( 1 \ I fu + ek 



X, 




(6) 
i B a J \ g u 



From the characteristic equation {A a — f u — eA a }{X a — g v — aeA a } — f v g u = 0, a pair of conjugate 
growth rates are obtained for each Laplacian mode as A a = (l/2){/ n + g v + (1 + a)eA a ± [4f v g u + 
(fu — 9v + (1 — &)£A a ) J 1 ' 2 }. Only the upper branch can become positive and it is always chosen 
as \ a in our analysis. From the condition that X a touches the horizontal axis at its maximum, 
the critical value of a is determined as a c = {f u g v - 2f v g u + 2[f v g u (f v g u - f u gv)] 1/2 }/fu an d the 



critical Laplacian eigenvalue is determined for a given e as A c = {(f u — g v )a c — \/\f v \guO'c(o'c + 
l)}/{ea c (a c — 1)}. The critical eigenvector in the (u,v) plane is given by (1,B C ) where B c = 
{-fu + gv + (cr c -^)eA c + [4f v g u + (f u — g v — (a c - l)(eA c )) I 1 / 2 }/(2/A These expressions coincide 



211]). if we replace there A by 



with the respective expressions for the continuous media (see, e.g., 
— q 2 , where q is the wavenumber of the plane wave mode. 
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FIG. 1: Linear stability analysis. Linear growth rates A a of the Laplacian modes a = 1,---,N for the 
Mimura-Murray model on a scale-free network (N — 200 nodes and mean degree (k) — 10) for three values 
of the diffusional mobility e and for the critical ratio of diffusion constants a — 15.5 ~ a c . Three curves 
corresponding e = 0.425, 0.165, and 0.060 are plotted as functions of the Laplacian eigenvalues A a . For 
comparison, curves with a — 15.0 and a = 16.0 are also drawn for e = 0.060. Critical modes are indicated 
for each value of e. The critical modes and the corresponding Laplacian eigenvalues are a c = 15, A c = —3.62 
for e = 0.425, a c = 135, A c = -9.32 for e = 0.165, and a c = 190, A c = -25.3 for e = 0.060. 



15 



(a) 0.6 

0.4 

.„ 0.2 

-©■ 




• 
« 

* 


0.0 1 


> r, o „ °W, ° 8w&*W#fteiiiMMf 


\. Q ° *P BI W$9P!S! 
__ ^ktP 3 ° 


-0.2 

n a 


1 . .1 1 . 



60 



40 



20 



2 3 

In ; 




(c) 



o 



o 



CD CP 



o 
oo 



o 



o 



o°co° O o 



o 



c° o 



o 



o°o^ OO o o ° 



o 



o o o o ' o o o o ooo 

%°o q.8 # oV.Oo 8 oo 

o r> • »o o o 

O Oq! A*»^° °^0 



o 



o 



o o 
o 



o 



o o» 



OO 



o 



o 



oo 



Oc9 



o 



o nn o "o c9 o o o o o 
n o o oogoOo8o Q 

O <O O ° ° n O 

°° o °oo° 
o °o ° °o 



O o 



o 



(d) 



•• 



o°«o° t o • • 
o oo o ° o 0(D0 . 

o - o o ooo oo o 0° 

dT o y °o o o too 
oo °o° o^o " 

o % 0o ^ 8 o°°o *° 00 o o 8 o -o 

- ■ O Q or pO O O q 

oo o o°° o 

'rvPn O C? O O „ • 

oo°ooSoOo8o # 



• - J o 
o • o 
o 



o 



- ft 

o o o 



,o o 



o 



oo co o 



o° ° 

,0° oo 



• o o oo u- - 

o . •o*o 



•• 



FIG. 2: Critical Turing modes for a scale-free network of size N = 200 and mean degree (k) = 10. (a,b): 
Critical eigenvectors (a) a c — 190 and (b) a c = 15 plotted against the node index i. Node degrees ki are 
shown by green stepwise curves. Node indices {i} are sorted according to their degrees {ki}. (c,d): The 
same critical eigenvectors (c) a c = 190 and (d) a c — 15 displayed graphically on the network. 
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FIG. 3: Localization of Laplacian eigenvectors in scale-free networks. The network size and mean degree 
are (a) N = 200, (k) = 10 and (b) N — 1000, (k) = 20. All nodes are divided into groups with equal 
degrees. For each group, the number of "differentiated" nodes with <pf' > 0.1 or <pf < —0.1 for each 
eigenvector a is counted. Then the fraction of such nodes in each group for each Laplacian eigenvector a 
is determined. Thus, these diagrams show density distributions of differentiated nodes for the entire set of 
Laplacian eigenvectors. 
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FIG. 4: Nonlinear evolution and a stationary Turing pattern of the Mimura-Murray model on a scale-free 
network at e = 0.12 and a = 15.6. The network size and mean degree are N = 1000 and (k) — 20. (a) The 
critical mode (the Laplacian eigenvector with a c = 422), (b) the activator pattern at the early evolution 
stage (t — 200), and (c) the stationary activator pattern at the late stage (t = 1500). Nodes are ordered 
according to their degrees; with (d) showing the dependence of the degree on the node index. 
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FIG. 5: Hysteresis and multistability. (a) Amplitude A of the Turing pattern vs. the diffusion ratio er; 
variation directions of a are indicated by arrows. The inset shows the blowup near R. (b) Stationary Turing 
patterns at the parameter points P (er — 17.0), Q (a — 13.5), and R (a = 12.8). 
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FIG. 6: Developed Turing patterns and mean-field bifurcation diagrams. Stationary Turing patterns com- 
pared with the bifurcation diagrams of the activator-inhibitor system on a single node coupled to global 
mean fields. The parameters are e = 0.12 and (a,b) a = 15.6, (c,d) a — 30. Blue curves (dots) indicate stable 
branches and light-blue curves (dots) correspond to the unstable branches. Crosses show the computed Tur- 
ing patterns. The global mean fields are (H (u \H (v ^ = (4.95, 9.97) for a = 15.6 and (#(*),#(*>) = (4.8, 9.9) 
for a = 30. 
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